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Abstract — Snow accumulation over an ice sheet is the sole mass 
input, making it a primary measurement for understanding the 
past, present, and future mass balance. Near-surface frequency- 
modulated continuous-wave (FMCW) radars image isochronous 
firn layers recording accumulation histories. The Semiautomated 
Multilayer Picking Algorithm (SAMPA) was designed and de- 
veloped to trace annual accumulation layers in polar firn from 
both airborne and ground-based radars. The SAMPA algorithm 
is based on the Radon transform (RT) computed by blocks and 
angular orientations over a radar echogram. For each echogram’s 
block, the RT maps firn segmented-layer features into peaks, 
which are picked using amplitude and width threshold parameters 
of peaks. A backward RT is then computed for each corresponding 
block, mapping the peaks back into picked segmented-layers. The 
segmented layers are then connected and smoothed to achieve a 
final layer pick across the echogram. Once input parameters are 
trained, SAMPA operates autonomously and can process hun- 
dreds of kilometers of radar data picking more than 40 layers. 
SAMPA final pick results and layer numbering still require a 
cursory manual adjustment to correct noncontinuous picks, which 
are likely not annual, and to correct for inconsistency in layer 
numbering. Despite the manual effort to train and check SAMPA 
results, it is an efficient tool for picking multiple accumulation 
layers in polar firn, reducing time over manual digitizing efforts. 
The trackability of good detected layers is greater than 90%. 

Index Terms — Antarctic ice sheet, image transforms, layers’ 
trackability, radar echo sounding, Radon transform (RT). 

I. Introduction 

R ADAR echo sounding of the ice sheets is and has been 
an active area of glaciological research, e.g., [1] and [2]. 
When a radar wave penetrates an ice sheet, the physical and 
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chemical properties of the firn reflect, refract, and attenuate 
the wave, creating an image of the internal structure (e.g., [3], 
vol II; [2]). Radar echograms have been used extensively to map 
the underlying bedrock of ice sheets, to determine past ice- 
flow histories and to determine accumulation rates, e.g., [4]- 
[8]. With the start of NASA’s Operation IceBridge in 2009, 
an unprecedented amount of new depth- sounding and near- 
surface radar data have been collected, and there is great need 
for automated and/or semiautomated algorithms for detecting 
and mapping isochronal layers in these data (e.g., [9] and 
[10]). There is particular need for automated methods that can 
resolve annual accumulation rates from the high-resolution and 
frequency-modulated continuous -wave (FMCW) radars that 
image the near-surface (« top 20-80 m) firn across the ice 
sheets (e.g., [11], [12]). Snow accumulation is the sole input 
to ice- sheet mass balance and is thus a necessary input for 
monitoring ice-sheet contributions to sea level rise. Reanalysis 
and regional atmospheric models are most commonly used 
to investigate changes in spatial and temporal variation in 
accumulation across ice sheets but are largely unchecked by 
measurements (e.g., [13]— [18]). Spatially and temporally exten- 
sive radar-derived measurements of accumulation are extremely 
important due to the paucity of accumulation measurements 
from ice cores and snow pits and the lack of a reliable annually 
resolved satellite retrieval of accumulation (e.g., [19]— [21]). 

Here, we present a first step toward retrieving spatially 
and temporally extensive maps of accumulation from radar 
measurements by introducing the Semiautomated Multilayer 
Picking Algorithm (SAMPA). SAMPA uses image-processing 
techniques to trace and extract firn layers in near- surface radar 
echograms. Layer extraction results are shown, along with 
a comparison of extracted layers to an ice core site in West 
Antarctica, showing the algorithm’s ability to detect annual 
snow accumulation. Although not presented here, SAMPA’ s 
layer extracting utility could be used to investigate other 
scientific questions involving near-surface layering, including 
layer deformation from wind or layer compaction due to 
densification. 

II. Background 

Tracking continuous isochronal internal layers over the ice 
sheets for hundreds to thousands of kilometers is difficult. Lay- 
ers can naturally bifurcate or become less pronounced, while 
noncontiguous layers may appear for short distances related to a 
local anomaly not associated with the broad climatic pattern or 
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event. Although difficult to extract, many approaches have been 
used to trace and track internal layers or ice surface and internal 
layers from radar echograms. Many studies have manually dig- 
itized layers (e.g., [7], [22], and [23]). This approach becomes 
prohibitively slow for long (thousands of kilometers) profile 
distances. The authors in [24]-[26] all used an interactive 
semiautomated means to insure accuracy and retrieve layers 
where discontinuities may occur. In [27], a method for the 
detection of near- surface Martian ice layers in orbital radar 
data was proposed. The method applies a series of filtering: 
first, a low-pass Gaussian filter enhances layers; next, a high- 
pass filter normalizes areas with high reflectivity with areas of 
low reflectivity. The authors then used a filter matched to the 
shape of the ice layers. The final layer detection is achieved by 
applying a threshold and morphological processing. Consider- 
ing our high-resolution data set approximates 4. 5 -cm vertical 
resolution x 20-cm along-track resolution, filtering steps com- 
bined with morphological processing are prohibited. Sime et al. 
[28] developed a fully automated processing method for picking 
internal layers and bed reflections using depth- sounding radar 
data. While the proposed method is useful for applications 
using depth- sounding radars, the resolution of the resulting 
extracted layers makes the application to high-resolution near- 
surface radars difficult. Specifically, the noise reduction step 
in the algorithm in [28] uses a horizontal and vertical moving 
average that degrades the original echograms’ horizontal and 
vertical resolutions, resulting in thicker extracted layers. A ma- 
jor challenge when using conventional edge detection processes 
is the location of the resulting edges (e.g., [29]). Commonly, 
the resulting edges’ location is fuzzy, which is an issue for 
applications trying to detect closely spaced annual layering. 
Very recently, Ferro and Bruzzone [30] have developed an 
automatic internal layer extraction method from radar sounder 
data sets. Their method uses several steps to carry the final layer 
extraction. In one particular step to highlight layer features be- 
fore extracting them, the original resolution of the image is lost 
and, therefore, the subsequent layer extraction location. This 
method lacks a regularization step; therefore, some extracted 
layer segments are not connected and, thus, their trackability 
making it less than ideal for our ice-sheet application. 

Continuous layers in near- surface radars, ranging from high 
to ultrahigh frequencies, have been shown in numerous studies 
to track isochronal layers over the ice sheets. The authors in 
[7] and [31] first showed that shallow 400-MHz radar data, 
which were gathered over the West Antarctic Ice Sheet (WAIS), 
contained isochronal layers with datable accuracy to ice cores 
of less than one year. In addition to high-/very high-frequency 
radars, ultra-/super-high-frequency (2-20 GHz), FMCW radars 
have been shown to map stratigraphic layers over ice sheets 
and ice caps at nearly annual rates (e.g., [32]-[34]). Addition- 
ally, Hawley et al. [35] first showed that an airborne radar at 
13.2 GHz (Ku-band) could determine annual layers in the 
dry- snow zone of the Greenland ice sheet when compared to 
annual density peaks. Radar-echo- sounding layers in the near 
surface arise from a natural seasonal change in snow density 
on an annual cycle. These studies prove the utility of near- 
surface radars for measuring accumulation, but all used manual 
time-consuming methods for picking isochronal layers. The 
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approach presented in our study to track and trace surface 
and internal layers from radar echograms is based on the 
Radon transform (RT) computed from the original vertical and 
horizontal resolutions of the echogram, thus preserving the 
high-resolution image of the firn microstructure. The advantage 
of using the RT is that the transformation is rooted from a 
summation along layer features and it has the capability of car- 
rying linear feature detection in a noisy surrounding with high 
location accuracy of the linear feature pixels, e.g., [36]. Unlike 
edge detection approaches that measure contrast between linear 
features and their surroundings, the RT maps linear features into 
a transformed domain, where thresholds are set for the peak’s 
amplitude and width detection. Our method detects segmented- 
layer features over small blocks, assuming that, over a short- 
enough distance, a curved layer feature may be approximated 
with a line segment, and thus validating the use of the RT. 
We have organized the remaining article as follows. Section III 
presents the RT applied to firn layer detection; the probability 
density function (PDF) of peaks in the RT domain is also 
analyzed. The different steps of the SAMPA algorithm are 
provided in Section IV. Results, analysis, and an example of 
application of SAMPA’s outputs are given in Section V. Finally, 
Section VI summarizes the SAMPA algorithm. 

III. RT Applied to Firn Layer Detection 

Radar-sounded isochronous layers within the ice sheets re- 
semble linear features carrying a variety of characteristics, 
including contrast inhomogeneity, discontinuity, natural bifur- 
cation, and speckle noise. A variety of approaches for extracting 
imaged linear features have been proposed. The Canny-Derich 
detector [37] behaves well on natural images where high con- 
trasts are measurable. Edges extracted using the Canny-Derich 
detector, however, suffer from a lack of location accuracy in 
certain scenarios (e.g., [29]). The ratio of local means has 
proven efficiency for amplitude radar imagery (e.g., [38]) and 
its applications, such as road detection (e.g., [39]). The RT (e.g., 
[40]) and its localized form (e.g., [41]) have been used for lines 
or line segment extraction in synthetic aperture radar images 
of rainforest areas (e.g., [36]). The summation inherent from 
the RT along linear features overcomes speckle effects in the 
resulting RT domain, where the detection of peaks is operated. 
While the conventional RT performs line detection even in 
noisy linear features, it cannot detect in the case of curved 
linear features naturally occurring within the isochronal ice- 
sheet layers. Parabolic and hyperbolic RTs have been used for 
certain specific applications (e.g., interpolation of seismic data 
[42], [43]), but do not replicate curved lines. Isochronous firn 
layers from ice sheets are randomly and continuously bending 
linear features, which are not ideal for a conventional RT. 

Here, to take advantage of properties of the conventional RT 
for layer detection, such as the ability to carry linear feature de- 
tection in a noisy image and the extraction of segmented linear 
features, we first process curved features within small enough 
blocks to be considered straight, and second, we perform 
integrations of curved features by summations using a short 
integration length. Thus, we compute the forward RT within the 
corresponding blocks, and the resulting RT domain is marked 
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Fig. 1. Synopsis of a line feature (LF) in the echogram space represented 
with (p, 6 , l) coordinates. The parameter l in (2) and (3) corresponds to 
the integration length LO (shown in the figure) along the linear feature, p 
and 0 are referenced from the center of the image. Elementary summations 
Si = Y, ^ p n are computed over LO, and point sums in the DRT domain 
correspond to N ■ Si, where N is the number of elementary summations 
performed over LO along the LF. 

by dark or bright spots, which correspond to dark or bright 
segmented-layer features representing the curved features. Spot 
peaks can be more easily detected in the transformed domain 
than the original curved features in the radar echogram space. 
The back projection of the peaks from the transformed domain 
allows the extraction of segmented-layer features. 

A. Mapping Multiple Firn Layers Into the RT Domain 

Radar echograms of ice sheets are distinguishable by strong 
radar signal returns from interfaces caused by annual density 
changes. The echograms image these annual layers. The use 
of the RT maps these layer features into bright spots in the 
RT domain. Thus, the RT domain is made by point sums. 
Bright spots in the RT domain are detected using the peak’s 
amplitude and width threshold parameters. That property of the 
RT makes it an efficient tool to highlight either dark or bright 
linear features even in a noisy and discontinuous surrounding. 

1 ) Forward RT (RT): The continuous RT 7£ c (p, 0) of a 2-D 
lattice I(x,y) with a compact support including the origin (the 
center of the lattice) is defined by [44] 

oo oo 

n c { P ,e) = II I (x,y)S(x cos 6 + ysinO — p) dxdy (1) 

— oo — oo 

where (x, y) corresponds to a given Cartesian position within 
the lattice space, p G (— oo, oo) is the distance from the lattice 
center to the line feature, 0 G [0, n) is the angle between the 
line feature and the horizontal axis (see Fig. 1 for illustration), 
and S is the Dirac distribution that converts the 2-D integral to 
the line integral along x cos 6 + y sin 0 = p. Hence, the RT is a 
set of projections along the angular directions 0 of all potential 
line features from the distance p to the center of the lattice. 

We consider that a point (x, y) of a given line feature making 
an angle 0 with the x-axis is specified by the following three 


real parameters: (/,p, 0), where l G (— 00 , 00 ) is a parameter 
along the line feature to integrate over it. Thus, the aforemen- 
tioned parameterization implies that, for a fixed direction 6 , the 
Cartesian coordinates (x, y) can be expressed in terms of local 
coordinates (/, p) on the line feature by 

( x = p cos 0 — Ism 0 
\ y = p sin 0 + / cos 6 

such that (1) can be rewritten after combining with (2) as 

00 

K c (p,6) = I I (p cos 6 — l sin 6, p sin 6 + l cos 6) dl (3) 

—00 

where dl denotes the line integral along the line feature 
x cos 0 + y sin 0 = p. 

Practical applications of the continuous RT deal with discrete 
domains such as images or echograms. Thus, we can write the 
discrete RT (DRT) 1Zd(p, 0) from (3) as 

00 

n d { P ,d) = I (p cos 6 — l sin/9, psmO + l cos 6) (4) 

i =— 00 

where the discretization step is considered to be one unit. 

2 ) Bright or Dark Spot Peak Detection in the DRT Domain: 
Fig. 2(a) and (b) depict a piece of an original echogram and the 
resulting DRT domain, respectively. Segmented-layer features 
in Fig. 2(a) correspond to peaks in the transformed domain in 
Fig. 2(b). The use of thresholds allows the detection of peaks’ 
amplitude and width, as shown in Fig. 2(c). Correspondingly, 
if a linear feature in the echogram domain is bright or dark, 
then the resulting spot in the DRT domain will be bright or 
dark and will represent a maximum or a minimum point sum 
value Si. Assuming that all point sums Si in the DRT domain 
are independent and identically distributed (supposing that the 
length of integration l = LO is small enough, i.e., LO ^ 1 m 
for relatively flat layers, for our application), the PDF PL(rn) 
of detecting a peak in the DRT domain (i.e., a segmented layer 
within blocks), if it exists, is given by [45] 

FL(m) 

TVdrt h(m)(f^h(m)drn) Nr>KT 1 , case of peak 

_ < (maximum) values 

Nnn T h(m)(l— J^°h(m)dm) JVdrt 1 , case of bottom 
, (minimum) values 

( 5 ) 

where h is the PDF of the distribution of the DRT domain, and 
TVdrt is the number of point sums in the DRT domain. 

The knowledge of the theoretical PDF of h is useful to 
analyze the behavior of the peak detection in the RT domain. 
Here, first, we have investigated the clutter, i.e., radar signal 
echo statistics modeling of the amplitude of Ku-band radar 
data of ice sheets. The goal was to derive a suitable PDF. The 
statistics modeling of the radar clutter distribution is important 
for algorithmic performance prediction, such as target detection 
and false-alarm probabilities [38], [39], [46]-[48]. Addition- 
ally, these statistics may be useful for Ku-band radar raw 
data processing consistency. To complete our investigation, we 
performed an analysis looking for the best theoretical fitted 
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Fig. 2. Line feature extraction using the DRT. (a) Piece of an original echogram exhibiting two line features corresponding to two firn layers hardly visible. We 
note discontinuous behavior and roughness of line features, (b) Resulting DRT domain. We note two main peaks, i.e., two point sums corresponding to the two 
line features, (c) Resulting threshold output of (b). (d) Result after the application of the backward DRT of (c). Here, l = L0 ~ 1 m (5 pixels) for relatively flat 
layers, PTV = 0.1, PWD = 4, and 0 E [— 7r/5, 7 t/5] . (e) Illustration of a vertical profile of PWD and PTV parameters for peak detection in the DRT domain. 
The greater the value of PWD, the more layers are detected; the smaller the value of PTV, the more layers are detected, (b) and (c) are in [x-axis = 0-axis and 
y-axis = p-axis], whereas (a) and (d) are in the original image reference, p and 0 parameters were sampled to the original image sizes in (a). 


model of the clutter Ku-band radar amplitude distribution of ice 
sheets. We compare the empirical PDF of Ku-band radar ampli- 
tude data with the following conventional and theoretical uni- 
modal distributions and their associated PDF for a given target 
(i.e., a pixel value p n > 0): log-normal, Weibull, gamma, and 
K distributions. We then used the conventional Kolmogorov- 
Smirnov test (KS-test) [49] and the Kullback-Leibler [50] 
distance (KLD) to carry out the goodness-of-fit theoretical 
distribution for the amplitude of Ku-band radar echograms. 
As a result, both the KS-test statistics and the KLD have 
substantiated that the theoretical log-normal PDF is the best fit 
for the amplitude of Ku-band radar echograms of ice sheets. 


Using a different data set, Ferro and Bruzzone [48] similarly 
studied the statistical distributions of radar sounder signals 
of ice sheets of the Shallow Radar (SHARAD) onboard the 
Mars Reconnaissance Orbiter of NASA. Their analysis sub- 
stantiated that the K distribution is suitable to model strong 
and weak layers. Here, the Ku-band radar data follow the the- 
oretical log-normal distribution £(p> 0 |/ij, 07 ) = (l / p<ji\/2it) 
exp(— ((ln(p) — pi) 2 /2crf)), where pi and 07 represent the 
mean and the standard deviation, respectively; therefore, a point 
sum Si of the DRT domain follows the sum of log-normal dis- 
tributions. Finding a close analytical expression of the sum of 
log-normal distributions, however, is an open problem; there is 
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Fig. 3. Theoretical probability of detecting peaks Pdp in the DRT domain 
versus the mean m (mu) and the standard deviation ai (sigma); with N = 10 
being the number of elementary summations over the length of integration 
LO = 5, A/drt = 50 x 50 the number of point sums in the DRT domain. Note 
that point sum values in the DRT domain are on the order of magnitude 10 9 ; 
hence, [ DRT ma x % = 10 9 , DRT max 2 = 10 10 ] is the range of point sums con- 
taining peaks in the DRT domain. On a given segmented layer within a block, 
the Pdp increases as 07 increases, and Pdp > 0.8, for 07 = 10 3 , for the range 
[DRT max 1 , DRT ma x 2 ] used for the simulation. This analysis substantiated 
that high signal return echoes (i.e., with a high variance of segmented-layer 
points) contribute to a high detectability of peaks in the DRT domain. 


no exact analytical expression. Using the following assumption: 
Within a block , pixels belonging to a given segmented layer 
have close intensity values po; therefore, a point sum in the 
DRT domain can be approximated by Si = N^2^ =1 p n ^N • 
LO • pq. It becomes much easier to assess the PDF I of a point 
sum Si (see more details in Appendix I) as follows: 


h{8i,N,L0,m,(r ,}ip) — 


N-L0 C (n '• J L0 >0 ^’ <J ') 


p<JiV2n 


exp - 


0 n ( aAo ) w) 


2°f 


( 6 ) 



Fig. 4. Description of SAMPA for ice-sheet radar echogram. Control param- 
eters are presented as follows: l = LO, the length of integration along line 
features, typically LO < 1 m (5 pixels) for relatively flat layers; PTV, the 
peak’s amplitude threshold value, which represents the height of peaks to be 
detected, typically PTV G [0.1, 0.5]; PWD, the peak’s width to be detected, 
corresponds to the width of peaks to be detected, typically PWD G [4,20]. 
B x size, By size, the block size dimensions, typically B x size x B y size < 
2.25 m x 20 m (50 x 50 pixels); 9, the angular orientation range, typically 
9 G [— 7r/2, — 7T / 2] ; and Ay « 22.5 cm (5 pixels) vertical distance within 
which layer points are considered to belong to the same layer feature. 


where N is the number of elementary summations Y^n = 1 Pn 
along a linear feature (see Fig. 1). 

From (5), the probability of detecting peaks Pdp , in the DRT 
domain, is then given by 


DRT max 2 

Pdp = J TL{m) 


dm 


DRT niax 1 

v, 


*DRT 




DRT, 


N-LO ) 

y/2 ) 


-erf 


aiV2 


AT-LO 


V2 


V2j 


(7) 


where TVdrt is the number of points [i.e., point sum S{] in the 
DRT domain, erf(x) represents the error function defined by 
erf(ar) = {2/^/n) f* exp (-t 2 ) dt , and [DRT maxl ,DRT max2 \ 
is the range of integration containing peak values in the 
DRT domain. Fig. 3 depicts the theoretical probability of 
detecting peaks Pdp in the DRT domain, for TVdrt = 50 x 


50 point sums, N = 10, L0 = 5, and the range [DLT max i = 
10 9 , LLI max 2 = 10 10 ]. On a given segmented layer within a 
block, the Pdp increases as cp increases, and Pdp > 0.8, for 
a 1 = 10 3 , for the range of point sums used for the simulation. 
This analysis substantiated that high signal return echoes (i.e., 
with high variance of segmented-layer points) contribute to a 
high detectability of peaks in the DRT domain, thus supporting 
the efficiency of using the RT to detect firn layer features from 
echograms. 


B. Firn Multilayer Feature Extraction From the DRT Domain 

After the detection of peaks in the DRT domain, due to the 
reversibility of the RT, a backward DRT is computed to get 
back only the expected segmented-layer features in the original 
space. There are a number of backward RT formulas. The 
following has been derived by [44]: 


7 r 00 

i{x - y)= PS j 


0 0 


{&RcDp)T&) 

p — (y cos 0 + x sin 0) 


dpde. (8) 






This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 



Fig. 5. Effect of variation of the length of integration Z = LO on extracted segmented-layers, B x size x B y size %10mx 2.25 m (50 x 50 pixels), PTV = 0.3, 
PWD = 20, and 9 E [— 7t/5, 7t/5], [a>axis = #-axis and y-axis = p-axis]. (a) Piece of an original amplitude of echogram: we note stronger return layers on the 
top and lower return layers on the bottom, (b) Log amplitude of echogram in (a). Note that the log amplitude better highlights layers comparing with (a) and was 
used, instead, to lower point sum values in the RT domain, (c) Extracted segmented-layers using LO « 1 m (5 pixels) for relatively flat layers and (d) extracted 
segmented-layers using L0~ 10 m (50 pixels); as l = L0 increases, we note three phenomena: smoothing of some layers, filling of some segmented layers, and 
depletion of some segmented layers. Smoothing and depletion cannot be corrected, whereas filling can be partially corrected using a regularization step. 


Conventional backward RT formulas are used to get the 
exact original distribution of I(x,y). Here, when focusing 
on peaks representing segmented-layer features in the DRT 
domain, after applying thresholds, we leave with few point 
sums [i.e., only peaks; see Fig. 2(c)] from which we would like 
to extract the original corresponding segmented-layer features. 
Therefore, the use of (8) is processing time consuming, and 
back projections to the original space suffice to extract the 
expected segmented-layer features using (4). 

Fig. 2(d) depicts the result of the application of the DRT, 
which is followed by the back projection of the two peaks in 
Fig. 2(c) to the original space of Fig. 2(a), allowing to get back 
the two corresponding segmented-layer features. 


IV. SAMPA for Ice-Sheet Radar Echograms 

Fig. 4 presents the different steps used in the SAMPA al- 
gorithm to track and trace isochronous fim layers from radar 
echograms. The input echogram is divided into blocks of same 
size, for instance, B x size x B y size , and each block is pro- 
cessed individually through the steps described in Fig. 4. Once 
each block is processed, the resulting extracted segmented- 
layers are projected back to the original block space. This 
approach was based on the assumption that a short-enough line 
segment, to even a curved line, can be considered straight and 
was also adopted to smooth the final layer pick. Note that the 
log amplitude of the imaged echogram better highlights layers, 
as shown in Fig. 5(b) compared to Fig. 5(a); moreover, the 
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DRT is a cumulative summation process. Applying the DRT 
on a log echogram (log amplitude) instead allows avoiding 
dealing with very large values (order of magnitude > 10 9 ) in 
the DRT domain. Thus, the log of the echogram is processed by 
following the steps. 

A. Step # 1 : Computation of the DRT 

The computation of the DRT requires specific parameters 
to be set. These parameters will vary based on the specific 
application and the quality of the radar echogram. The param- 
eters and their function in the SAMPA algorithm are described 
in the following. B x size and B y size represent the optimal 
horizontal and vertical dimensions of a block size, respectively; 
LO represents the optimal length of integration along the line 
features; and 0 represents the angular orientation range. The 
outputs of this step are DRT maps [see illustrations in Figs. 2(b) 
and 6(a)]. For the Ku-band radar echograms processed here, the 
optimal ranges for the three input parameters are the following. 

• The block size B x size x B y size is chosen within ~ 
[10 m x 2.25 m, 20 m x 4.5 m] ; using a much larger block 
size will result in smoothing the firn layers’ undulations 
into straight line- segments. Smaller block sizes will result 
in a long processing time and may increase false detections 
in the scenario where few radar echoes are available within 
blocks. 

• The choice of the integration length LO is also impor- 
tant. Smaller values of LO lead to fine integration along 
segmented-layer features for the computation of the DRT 
and the backward DRT, whereas greater values of LO 
lead to rough integrations. Smaller values of integration 
length, moreover, lead to a good reconstruction of the 
layer features, as shown in Fig. 5(c) and (d). A length of 
integration LO « 1 m for relatively flat layers is optimal 
for our application for the computation of DRT and the 
backward DRT of input echograms’ blocks. Increasing 
LO leads to the smoothing of some final picked layers, 
over filling of other segmented layers (i.e., possibility 
of connecting different spaced segmented-layer features 
by integration), and depletion of some segmented-layer 
features (i.e., case of presence of few layer echoes that 
can be integrated with non layer echoes). Smoothing and 
depletion cannot be corrected in any other step of the 
algorithm, but correcting the segmented layers is partially 
corrected in the regularization step. 

• Last is the choice of 6, which is the angular orientation 
range. 6 depends on the layers’ local slope. The flatter the 
layers, the smaller the angular orientation range. Imaged 
layers from echograms of our data set are relatively flat; 
the angular orientation range of [— 7 t/5,7t/ 5] works well. 
Using a much bigger range comparing to the layers’ local 
slope, however, will increase the DRT time processing, but 
may not improve the overall performances. 

B. Step #2: Peak Detection in the DRT Domain 

In Section III-A2, we demonstrated that high signal return 
echoes contribute to a high detectability of peaks in the DRT 


domain. High signal return echoes correspond to layer points in 
the echogram and peaks in the DRT domain. Therefore, we have 
high confidence that detecting peaks in the DRT domain leads 
to detecting corresponding layers in the echogram. Moreover, 
setting thresholds while detecting peaks in the DRT domain 
corresponds to fixing a certain level of probability of detection 
of segmented layers. Inputs for peak detection require threshold 
values to determine which maxima in the DRT domain should 
be picked as segmented layers. The inputs PTV, i.e., the 
peak’s amplitude threshold value, and PWD, i.e., the peak’s 
width to detect in the DRT domain, are needed for the detection 
of peaks, as illustrated in Fig. 2(e). The PTV describes the 
heights of the peaks, and the PWD describes the widths of 
the peaks. The outputs of this step are the peaks’ maps, as 
illustrated in Fig. 6(b) and (c). Greater values of PWD allows 
the detection of more peaks, whereas smaller values of PWD 
permits the detection of fewer peaks, as shown in Fig. 2(e). 
Between the two threshold parameters PTV and PWD, the latter 
may affect the result if not cautiously tuned. For our application 
using near- surface radar echograms with relatively flat layers, 
the following range is chosen by tuning the algorithm to most 
closely pick the visible layers in the echogram, i.e., PWD > 4. 
Moreover, the parameter PTV represents the percentage of the 
peak value to consider for the extraction of layer features from 
the DRT domain. As PTV increases, the number of peaks 
detected in the DRT domain decreases, leading to the extrac- 
tion of fewer layers. Smaller values of PTV not only allow 
the extraction of more layers but also tend to increase false 
detections. Layers’ signal return echoes are inhomogeneous, 
i.e., layer feature echoes do not have the same peak’s amplitude 
across the echogram. Thus, that creates an uncertainty while 
detecting peaks in the DRT domain and by, while extracting the 
corresponding segmented layers. That uncertainty leads to false 
detection. Moreover, the existence of uncertainty implies that 
the false detection rate cannot be null. A compromise has to be 
made based on the end-user goals to minimize false detections. 
The final results presented in this paper used PTV £ [0.1, 0.5]. 
Fig. 6(d) (PWD = 20 and PTV = 0.1) and Fig. 6(e) (PWD = 4 
and PTV = 0.1) depict the profile from the column where the 
peaks are located in Fig. 6(b) and (c), respectively, highlighting 
the effect of the variation of the two parameters. 

C. Step #3: Segmented-Layer Feature Extraction 

The detection of peaks from the DRT domain produces the 
segmented-layer features to be extracted. Thus, back projec- 
tions are computed for each peak leading to the extraction of 
segmented-layer features in the original space. The following 
inputs are set in this step: the block size dimensions B x size 
and B y size ; the length of integration LO; and the angular 
orientation range 0 used in Step #1. The output of this step is 
the map of segmented-layer features, as shown in Fig. 7(b). 

D. Step #4: Regularization of Segmented-Layer Features 

Fig. 7(a) depicts an original Ku-band radar echogram, and the 
resulting segmented-layer features are shown in Fig. 7(b). Some 
extracted segmented-layer features after applying the backward 
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Fig. 6. Effect of variation of PWD in the DRT domain for fixed PTV, LO « 1 m (5 pixels) for relatively flat layers, and 0 E [— 7r/5, 7t/5], [a>axis = 0-axis and 
y-axis = p-axis]: (a) DRT domain computed from a piece of an original amplitude echogram in Fig. 5(a). (b) PTV = 0.1 [i.e., > 10% of the maximum value 
in the DRT domain in (a)] and PWD = 20. (c) PTV = 0.1 and PWD = 4. (d) Profile of the column where the DRT maximum values are detected, PWD = 20 
and PTV = 0.1. (e) Profile for PWD = 4 and PTV = 0.1. For a fixed PTV, as PWD increases, the number of peak values detected in the DRT domain increases 
leading to the detection of more layers [note the number of bullet peaks in (d)]. Smaller values of PWD allow the detection of fewer layers [see in (c)]. 


DRT are noisy and disjunct, as shown in Fig. 7(b). That is due 
to detected nonphysical layers and to partly contrasted layers 
mostly resulting to anomalies. Visible firn layers on echograms 
fit into different categories, including Contrasted/Less Con- 
trasted Layers Across the Echogram (CLAE/LCLAE) and 
Partly Contrasted/Less Contrasted Layers Across the Echogram 
(PCLAE/PLCLAE). We do not consider compressed layers 


here, which may be hard visible. In some scientific appli- 
cations such as accumulation estimation, the trackability of 
most extracted layer features is required. Therefore, underlying 
segmented layers of layer features that have lost parts due 
to anomalies need to be connected. Thus, resulting extracted 
segmented-layers are treated based on their trackability across 
the echogram as follows: CLAE/LCLAE and PCLAE/PLCLAE 
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Fig. 7. Layer tracking results, (a) Original Ku-band radar data 03 echogram 0008. (b) Extracted segmented-layers. (c) Final result after the smoothing procedure, 
(d) Final result overlapping the original echogram in (a). B^size x B^size % lOmx 2.25 m (50 x 50 pixels); L0 « 1 m (5 pixels) for relatively flat layers; 
PTV = 0.1; PWD = 20; and 6 £ [— 7r/5, 7 t/5] . Upper layers are better extracted than lower layers due to lower signal returns in deeper layers. Thus, the 
reliability of the extraction of deeper layers decreases along with their trackability. 


lying within less than half the long track of the echogram are half the long track of the echogram are kept and then connected 
ignored due to their lack of trackability across the echogram, and smoothed. Hence, a regularization process that aims at con- 
CLAE/LCLAE and PCLAE/PLCLAE lying within more than necting and smoothing extracted segmented-layers is necessary 
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TABLE I 

All Parameters Used Along With Their Functions in the SAMPA Algorithm 


Parameters 

SAMPA Step set in 

Functions 

B x size , By size 

DRT and backward DRT 

Block size: used in the computation of the forward and the backward DRT by block 

L0 

DRT and backward DRT 

Length of integration: used in the computation of the forward and backward DRT 

e 

DRT and backward DRT 

Orientation angle range: used in the computation of the forward and backward DRT 

PTV 

Peak detection in DRT space 

peak threshold value: sets the minimum amplitude for a peak to be considered 

PWD 

Peak detection in DRT space 

peak’s width to detect: sets a range of amplitude values within which peaks are considered 

Ay 

Smoothing process initialization 

Vertical bins distance within to consider two segmented-lines belonging to the same layer 


to achieve the final layer tracking. Extracted segmented-layers 
are then connected if they are vertically within Ay pixels and 
are indexed to be the same feature. Tests have conducted to 
choose Ay « 22.5 cm for echograms of our data set. The use 
of lesser values lowers the layer detection probability, whereas 
greater values tend to connect return signal echoes from dif- 
ferent depth. The regularization step connects and smooths 
extracted layers, as depicted in Fig. 7(c). There are numerous 
smoothing algorithms. Here, we use a smoothing procedure 
based on the fully automated algorithm developed in [51]. The 
smoothing algorithm can work with evenly spaced data, deals 
with weighted data, can account for the occurrence of missing 
values, and self-carries the estimation of the smoothing pa- 
rameter. The smoothing algorithm estimates the final smoothed 
layer from the extracted segmented-layers using an iterative 
scheme algorithm (see Appendix II for details). In practice, 
less than ten iterations are sufficient to get an acceptable 
smoothing result. 

All parameters used by the SAMPA algorithm and their 
functions are described in Table I. Although the DRT is a ge- 
ometric integral, the resulting DRT domain carries radiometric 
information from the original echogram due to the existence of 
the inverse RT. Hence, the radiometry information is taken into 
account while detecting peaks in the DRT domain. However, 
peak detection in the DRT domain carries uncertainty inherent 
from the radiometric layer features’ uncertainty of the original 
echogram, uncertainty that results to false or missing detection. 

E. Step #5: Averaging of Layers’ Picks Along Track 

Once the picked vertical layers are established by SAMPA, 
the picks are summed along the y - axis for some given distance 
along the x-axis, in our case km along track. This results in 
a single layer depth pick across the specific x distance. SAMPA 
uses the higher resolution radar data for all of the layer picking, 
but this final step reduces the resolution specifically for the final 
physical application. This averaging step could be turned on or 
off depending on the scientific application. 

V. Results and Analysis 
A. Data Description 

We tested the SAMPA algorithm on two Ku-band radar 
echogram data sets, i.e., one airborne and one ground-based: 
The Ku-band Radar was built by the Center for Remote 
Sensing of Ice Sheets (CReSIS) at the University of Kansas 


(https://data.cresis.ku.edu/#KBRA), and it has a typical sweep 
from 13 to 17 GHz. This radar images the near- surface firn 
to depths of over 20 m, in dry firn conditions. The technical 
specifications of the instrument are given in [52], and the data 
processing detailed is outlined at (http://nsidc.org/data/docs/ 
daac/icebridge/irkublb/). The primary purpose of this radar 
is high-precision surface elevation measurements over polar 
ice sheets, although the instrument can be also used to map 
annual layers with a high vertical resolution of « 4.5 cm 
and along-track resolution of ~ 20 cm. The radar illuminates 
individual spots, with a diameter size of « 2 m, that are stacked 
to improve the along-track resolution. Hence, the echogram 
pixel size approximates 4.5 cm x 20 cm. The original data 
domain is the amplitude. To lower DRT domain values, the data 
are then converted into decibels before processing. Here, we 
use echograms collected during two different field campaigns 
that used the same Ku-band radar. SAMPA was run on an 
airborne data set, i.e., the NASA’s IceBridge Ku-Band Radar 
LIB Geolocated Radar Echo Strength Profiles data, collected 
during one flight over West Antarctica by the National Snow 
and Ice Data Center [53]. We also tested SAMPA on a ground- 
based data set collected during the Satellite Era Accumulation 
Traverse (SEAT) in 2010 [54]. Here, we show only the results 
from the ground-based data, which are similar to the airborne 
results. The ground-based data were used because it has a 
higher along-track resolution of 20 cm versus 5.6 m for an 
airborne data set. 

B. Results 

The SAMPA algorithm was applied to a radar echogram 
taken in West Antarctic. This region is known for relatively 
flat firn layers. The algorithm was run and the parameters 
trained. For our application using near-surface radar echograms 
with relatively flat layers, the following values were chosen by 
tuning the algorithm to closely pick the visible layers in the 
echogram: the block size B x size x B y size ^10mx2.25m 
(50 x 50 pixels); the integration length LO « 1 m (5 pixels); 
the peak’s width to detect PWD = 20; the peak’s amplitude 
threshold value PTV = 0.1; the angular orientation range 6 = 
[— 7r/5, 7t/5] ; and Ay ^22.5 cm (5 pixels vertically). 

For the results reported in this paper, we discarded CLAE/ 
LCLAE and PCLAE/PLCLAE lying within less than < 200 m 
of the echogram, due to their lack of trackability across wide. 
However, CLAE/LCLAE and PCLAE/PLCLAE lying within > 
200 m of the echogram are retained for the regularization step. 
Fig. 7(c) depicts the final result of extracted layers from the 
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Fig. 8. Additional layer tracking results: (al) original Ku-band radar data 03 echogram 0020; (a2) original Ku-band radar data 02 echogram 0068; (bl) final result 
data 03 echogram 0020 after the smoothing procedure of al); (b2) final result data 02 echogram 0068 after the smoothing procedure of a2); B x size x B y size « 
10 m x 2.25 m (50 x 50 pixels), l = L0 « 1 m (5 pixels) for relatively flat layers, PTV = 0.1, PWD = 20, and 6 E [— 7t/5, 7t/5] . Bottom layers in (a2) are 
much less pronounced, therefore, the resulting extraction in (b2). 


TABLE II 

Summary of All Parameter Settings Used by the SAMPA Algorithm. Recommended Ranges and Values Were 
Used for Results Presented in the Study. B x size x B y size: Block Size’s Dimensions, L0: Integration Length, 
and 6: Angular Orientation Are Used to Compute Forward and Backward DRT. PTV: Peak’s Amplitude 
Threshold Value, PWD: Peak’s Width to Detect, and Ay: Vertical Shift Within Which Two Segmented-Layer 
Features Are Considered to Belong to the Same Layer Are Used in the Detection Step. Here, B x size, 
B y size, L0 Are Expressed in Number of Pixels, and a Pixel’s Dimension « 4.5 cm x 20 cm 


/ 

Forward DRT 

Spot Peak detection 

Backward DRT 

Parameters 

B x size , B y size 

L0 

0 

PTV 

PWD 

Ay 

B x size , B y size 

lo e 

Range 

[50, 100] 

[1,5] 

[-tt/2, tt/2] 

[0.1 0.5] 

[4,100] 

[5,20] 

[50,100] 

[1,5] [-tt/2, tt/2] 

Value used 

{50} 

{5} 

[-71-/5, 7t/5] 

{0.1} 

{20} 

{5} 

{50} 

{5} [-tt/5, tt/5] 


Ku-band radar echogram in Fig. 7(a) after the regularization 
step. Most picked layers are continuous and smoothed. There 
is no better way of assessing the localization accuracy of an 
extracted layer feature than lining it up above the original 
echogram. Hence, we superimposed extracted layers with the 


original echogram in Fig. 7(d) for a visual check of location 
accuracy of picked layers. Additional layer extraction results 
are shown in Fig. 8(bl) and (b2). Table II summarizes all 
parameters and recommended ranges and their exact values set 
for the final results presented in this paper. 
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TABLE III 

SAMPA Layer Tracking Performance Evaluation From the Echograms in Fig. 7(a) and Fig. 8(al) and (a2). These Performances Are 
Indicative and May Lightly Vary Based on the Quality of the Echogram. #X, X = Performance Number and the Corresponding 
Percentage; #V, V = Visually Detected Layers by an Expert; #G, G = Good Detected Layers, i.e., Layers That Are Visually 
Detected by an Expert and Are Automatically Detected by SAMPA, Including Partially Detected Layers; #F, F = False 
Detected Layers, i.e., Layers Detected by SAMPA but Are Not Visually Detected by an Expert; #M, M = Missing Detected 
Layers, i.e., Layers That Are Visually Detected by an Expert but Are Not Automatically Detected by SAMPA. Columns Layer 
SAMPA: Good and Layer SAMPA: False Represent Tracking Rate and False Detection Rate Relatively to Visually Detected 
Layers by an Expert, Respectively. We Note That SAMPA Provides Better Detection Rate and Lower False Detection 
Probability for Upper Layers Than for Deeper Layers. Global Performances Are Also Provided for Each Echogram 


UPPER Layer Tracking LOWER Layer Tracking 


Data 03 Echogram 0008 - Global; V =60/60 (100%), G=82%>F=17%,M=18% 


Visual: #V 

SAMPA; #G 

SAMPA: #F 

SAMPA: #M 

Visual: #V 

SAMPA: #G 

SAMPA: #F 

SAMPA: #M 

30 

30 

4 

0 

30 

19 

6 

11 

#v/#v 

#G/#V 

#F/#V 

mm 

#V/#V 

#G/#V 

#F7#V 

#M/#V 

30/ JO (100%) 

30/30(100%) 

4/30 (13.4%) 

0/30 (0%) 

30/30 (100%) 

19/30 (63%) 

6/30 (20%) 

1 1/30 (37%) 


Data 03 Echogram 0020 ■ 

Global: V=57/57 (100%), G=84%, F=17.5%, M 

=16% 


Visual: #V 

SAMPA: #G 

SAMPA: #F 

SAMPA: #M 

Visual: #V 

SAMPA: #G 

SAMPA: #F 

SAMPA: #M 

27 

26 

3 

1 

30 

22 

7 

8 

#v/#v 

#G/#V 

#F/#V 

#M/#V 

#v/#v 

#G/#V 

#F/#V 

#M/#V 

27/27 (100%) 

26/27 (96%) 

3/27(11%) 

1/27 (4%) 

30/30(100%) 

22/30 (73%) 

7/30(23%) 

8/30 (27%) 


Data 02 Echogram 0068 

- Global: V=60/60 ( 100%),C=87%, F=28%, M= 

13% 


Visual: #V 

SAMPA: #G 

SAMPA; #F 

SAMPA: #M 

Visual: #V 

SAMPA: #C 

SAMPA: #F 

SAMPA: #M 

39 

34 

12 

5 

21 

18 

5 

3 

#V/#V 

#G/#V 

#F/#V 

#M/#V 

#V/#V 

#G/#V 

#F/#V 

#M/#V 

39/39 (100%) 

34/39 (87%) 

12/39 (31%) 

5/39(13%) 

21/21 (100%) 

18/21 (86%) 

5/21 (24%) 

3/21 (14%) 


C. SAMPA Layer Tracing Performance Evaluation 

We assessed the layer tracking capability of the SAMPA 
algorithm using the results shown in Fig. 7(c) and Fig. 8(bl) 
and (b2). From original echograms, upper layers are more 
contrasted and visible than lower layers due to lower signal 
returns from deeper layers. We, therefore, split each echogram 
into upper and lower portions from the middle, i.e., approxi- 
mately at 20 m depth, where the layers’ contrast starts dropping 
for evaluation in terms of percentages of the following: good 
detected layers (G), i.e., layers that are visually detected by 
an expert and are automatically detected by SAMPA, including 
partially detected layers; false detected layers (F), i.e., layers 
detected by SAMPA, but that are not visually detected by an 
expert; and missing layers (M), i.e., layers that are visually 
detected by an expert, but are not automatically detected by 
SAMPA, as highlighted in Table III. Thus, SAMPA provides 
a better layer detection probability and a lower false detection 
probability for upper layers than for deeper layers. In Table III, 
the global good detected layers’ rate G is greater than 80%, 
which is compatible with the theoretical PDF substantiated 
in Fig. 4. Our quantitative evaluation showed that, generally, 
SAMPA has good performance while extracting upper lay- 
ers, and thus, upper layers were much easier to track within 
echograms, which was expected due to the attenuation of the 
radar signal in deeper layers and their low contrast manifested 
by low variance, as proven by the theoretical PDF in Fig. 4. 
This quantitative analysis is indicative, since the final extracted 
layers have to be validated by the user according to its applica- 
tion. There is a possibility of connecting fake segmented layers 
during the running of the smoothing procedure. That may occur 


particularly for bottom layers where the contrast has been lost 
due to anomalies, including layer compression. That is a typical 
case where the SAMPA algorithm fails by vertically connect- 
ing nearby segmented layers. Fake layer connection, however, 
occurs when there are missing segmented layers within verti- 
cally 2 A y across the echogram. Hence, considering individual 
layers, fake segmented-layer connection are manageable by the 
end user according to the final scientific application. 

1 ) Parameters ’ Sensitivity and Performance Behavior: The 
detection of peaks in the DRT domain is a crucial step in the 
SAMPA algorithm. It encompasses the peak’s amplitude and 
width detection. The reliability of the resulting maximum spot 
map, from which the backward DRT is computed, depends on 
that step to produce extracted layers. Thus, the detection of 
peaks uses two threshold parameters PTV and P WD to high- 
light peaks in the DRT domain. Consequently, false detections 
in the final result mainly come from the tuning of those two 
parameters. Therefore, it is necessary to analyze the sensitivity 
of PVT and PWD on performance measures G, F, and M. 
In particular, we would like the percentage of false detected 
layers F to be low or at least constant for a given set of PVT 
and PWD values. Plots of percentages of good detected layers, 
of false detected layers, and of missing layers versus PWD 
are shown in Fig. 9(al)-(cl), respectively, for the echogram 
in Fig. 7(a), for different values of PVT. In addition, similar 
plots versus PVT are shown in Fig. 9(a2)-(c2), respectively, for 
the same echogram, for different values of PWD. Percentages 
of good detected layers in Fig. 9(al) increase for different 
values of PVT as PWD increases, and they stay constant in 
Fig. 9(a2) as PVT increases for different values of PWD. 
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Fig. 9. Parameters’ sensitivity and performance behavior. Echogram in Fig. 7(a) was used to analyze the sensitivity of parameters PVT and PWD for the 
performance behavior. For different values of PVT, as PWD increases: (al) percentages of good detected layers (G); (bl) percentages of false detected layers (F); 
(cl) percentages of missing detected layers (M). As PVT increases and for different values of PWD: (a2) percentages of good detected layers (G); (b2) percentages 
of false detected layers (W); (c2) percentages of missing detected layers (M). In (al), percentages of good detected layers increase for different values of PVT as 
PWD increases, and in (a2), they stay constant as PVT increases for different values of PWD. In (bl) and (b2), percentages of false detected layers stay roughly 
within [15%, 27%], for whichever PVT or PWD increases for different values of either PWD or PVT. In (cl), percentages of missing detected layers decrease for 
different values of PVT as PWD increases, and in (c2), they are constant when PVT increases for different values of PWD. This analysis substantiates that the 
SAMPA algorithm has a quasi-constant false detection rate relatively to visually detected layers by an expert, below 20% in general, and that overall performances 
are good for the following settings: PWD = 20 and PVT E [0.1, 0.5], for which the detection rate is around 80%. 


Percentages of false detected layers stay roughly within [15%, 
27%] for whichever PVT or PWD increases for different values 
of either PWD or PVT in Fig. 9(bl) and (b2). Percentages 
of missing detected layers decrease in Fig. 9(cl) for different 
values of PVT as PWD increases, and they are constant in 
Fig. 9(c2) as PVT increases for different values of PWD. 
This analysis substantiates that the SAMPA algorithm has a 
quasi-constant false detection rate, below 20% in general, and 
that overall performances are good for the following settings: 
PWD = 20 and PVT E [0.1, 0.5], for which the detection rate is 
around 80%. 

2 ) Trackability of Good Extracted Layers: One of the ulti- 
mate purposes of a layer extraction is the trackability of the 
resulting layers, i.e., being able to follow each individual layer 
across wide. In practice, two scenarios occur. 

• In an ideal scenario, extracted layers do not cross, and 
we have a perfect extraction, as shown in Fig. 10(al). 
If we weigh each extracted layer point to one unit, the 
sum of extracted layer points for each column is identical 
to the number of layers, as illustrated in Fig. 10(bl). 
The resulting characteristic appears horizontally, meaning 
perfect trackability of extracted layers. 

• In a common scenario, extracted layers may cross, may 
carry false detection, or may have missing detection, or 
local anomalies, as shown in Fig. 10(a2). In this case, the 
sum of extracted layer points for each column is different 


to the number of layers, as illustrated in Fig. 10(b2). 
The resulting characteristic displays downward extrema, 
meaning lack of perfect trackability of extracted layers. 

We have the characteristic of layer trackability (COT) as 
COT = Ececolumn numbers c), where L corresponds to the 
array of extracted layer points, where each layer point weighs 
one unit. Here, the way to interpret COT is about how well good 
detected layers are trackable ?, knowing that there exists few 
false detections. 

Understanding COT interpretation is simple, but defining a 
COT index iCOT , i.e., an index of goodness of good detected 
layers’ trackability, that tells whether extracted layers are well 
trackable or not across wide is a hard task. Indeed, such a 
COT index interpretation must take into account additional 
information including a global criterion, such as good and false 
detection rates, to avoid fully trackable layers but maybe made 
only of false layers; but also a local criterion such as measure 
of trackability per column wide across the echogram given 
the possibility that layers may cross, or that there are local 
anomalies leading to partly detected layers that are independent 
of false detection rate. 

Hence, we define here the low bound of COT index 

iCOTuin by 

iCOT Min = ^ mcOT e [0, 100%] (9) 

Max cot 
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Fig. 10. COT: Here, only good detected layers were considered, (al) Ideal scenario, a piece of extracted layers with perfect detection. (a2) Common scenario, 
a piece of extracted layers with crossing layers, (bl) COT from (al) is a horizontal line, meaning perfect trackability, the characteristic of layer tracking index 
iCOTMin = 100%, iCOT Avg = 100%, and the validation criterion of iCOT, VCqot = 100%. Here, ideal case = 0. (b2) COT from (a2) appears with 
downward extrema (upside-down pyramid), meaning lack of perfect trackability, VCqot = 87.5%, iCOT M [ n = 85.7%, iCOT Avg = 97.62%; here, = 
3/0. (c) COT in Fig. 7(c) (VCcot = 83%, iCOT M in = 92%, iCOT Avg = 98.95%, #F = 10). (d) COT in Fig. 8(a2) (VC G o t = 83%, iCOT M in = 
92%, iCOT Avg = 98.85%, #F = 10). (e) COT in Fig. 8(b2), (VC c ot = 75%, iCOT M in = 93%, iCOT Avg = 98.98%, #F = 17). Values of iCOT < 
100% are due to crossing, or missing good detected layers, and we have a relatively good trackability index. More than 90% of good detected layers are trackable, 
as highlighted in Fig. 11. 


where MincoT and Mclxqot represent the minimum and the 
maximum of the COT’s values, respectively. We also define the 
average of COT index iCOT Avg by 

iC OT Avg = ^ gcOT G [0, 100%] (10) 

Max COT 

where AvgcoT corresponds to the average of the COT’s values 
column wide. COT here is computed only for L(:,c) repre- 
senting good detected layers. Therefore, any decrease of either 
iCOTMin or iCOT Avg comes from artifacts, including crossing 
or missing parts of good detected layers. Either (9) or (10) 
is dependent of the existence of the aforementioned artifacts. 
Indeed if just a column of good detected layers in L(:, c) is full 


of zeroes, iCOTyiin = 0, saying that we face bad trackability. 
However, iCOT Avg will be different to 0 as long as there are 
good detected layers. In the same spirit, one crossing of good 
detected layers decreases iCOTyiin - However, in the meantime, 
a good detected layer rate may be great, and in an operational 
stand point, layers with few missing detected points on columns 
or few crossings can still be trackable. Thus, iCOTyiin gives 
the worst scenario, whereas iCOT Avg gives average behavior 
of iCOT. 

Therefore, either (9) or (10) is not meaningful alone for all 
detection result scenarios. They are more meaningful in the case 
of a perfect trackability of good detected layers or for an indi- 
vidual COT’s column. However, an individual COT’s columns 
by themselves are less significant for a global interpretation 
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TABLE IV 

Behaviors of the Validation Criterion of COT VCcot, COT 
Minimum Index iCOT M i n > COT Average Index iCOT Av g , 
and False Detection Rate F for the Three Echograms 
Processed in the Current Study. zCOT M i n Values 
Correspond to Worst Scenarios, Whereas 
iCOT Avg Values Represent Average Behaviors 


/ 

Data 03 Echogram 0008 

Data 03 Echogram 0020 

Data 02 Echogram 0068 

V CcOT 

83 % 

83 % 

75 % 

iCOT]\/[ in 

92 % 

92 % 

93 % 

iCOT^yg 

98 . 95 % 

98 . 85 % 

98 . 98 % 

F 

10(17%) 

10(17.5%) 

17(28%) 


of goodness of trackability of good detected layers. We need 
an additional criterion that captures layer detection global per- 
formance first and then uses it to validate the interpretation of 
either (9) or (10). If we define a validation criterion (i.e., layer 
detection global performance) V Cqot as 

VCcot = #F e [o, 100%] (1 1) 

where #G and represent the number of good detected 
layers and the number of false detected layers, respectively. In 
an ideal scenario, given in Fig. 10(bl), = 0, VCcot = 

100%, MincoT = M axcoT, iCOTMin = 100%, and 
iCOT Avg = 100%; otherwise, as shown in Fig. 10(b2), 
i^F = 3 % 0, VCcot = 87.5%, MincoT < Max cot, 
iCOTuin = 85.7%, and iCOT Avg = 97.62%; both iCOT U in 
and iCOT Avg are less than 100%, meaning lack of 
perfect trackability. The higher the values of VCcot and 
iCOT M in,Avg> the more meaningful iCOT M in,Avg will be. 
iC OTjviin, Avg denotes either iCOTyiin or iCOT Avg • For a 
general validation purpose and based on performances given 
in Table III, if we set up VCcot > 70%, which means F 
would be small enough (< 30%) and G would be high enough 
(> 70%), then the use of 7COTMin,Avg as a trackability index 
in conjunction with VCcot would be meaningful. COTs for 
good detected layer results in Figs. 7(c), 8(bl), and (b2) are 
depicted in Fig. 10(c)-(e), respectively. 

Based on our assumption (only good detected layers 
with few short missing parts and few crossings are con- 
sidered for a meaningful iCOTMin,Avg), the correspond- 
ing (VCcot CCOTuin, iCOT Avg , F) values are given in 
Table IV. Values of VCcot, iCOT M in, and iCOT Avg are less 
than 100%, due to some crossing or missing parts of good 
detected layer points. Some layers may partly split by the end of 
the echogram, whereas others may partly join by the beginning 
of the echogram. That affects VCcot, iCOTMin, iCOT Avg , 
and F behaviors. 

COT plots resemble upside-down pyramids [see 
Fig. 10(b2)-(e)] showing different levels of MincoT values. 
MincoT values depend on the presence of artifacts more than 
Max cot- If Max cot values are considered #level = 0, the 
closer the #level to Max cot values, the higher the iCOTMm- 
The further the #level to Max cot values, the lower the 
iCOTMin, as the trackability of good detected layers becomes 
more sensitive to artifacts, as illustrated in Fig. 11. Thus, for 
the three examples presented in this paper, more than 90% 
of good extracted layers are trackable across wide, with a 



Fig. 11. iCOTmiu versus #levels of COT pyramid in Fig. 10(c) (+), 
Fig. 10(d) (o), and Fig. 10(e) (*). If Max cot values are considered #level = 
0, the closer the #level to Max cot values, the higher the iCOT M [ n . The fur- 
ther the #level to Max cot values, the lower the iCOTMin, as the trackability 
of good detected layers becomes more sensitive to artifacts, including crossing 
and missing parts of good detected layers. 

good detection rate VCcot greater than 70%. Layers in our 
area of interest are roughly flat; hence, our approximation of 
applying the DRT over small blocks holds. As we highlighted 
it in Section III- A, the DRT maps linear feature radiometry 
into peaks in the DRT domain where the detection operates. 
Remote sensing information carries uncertainty when it comes 
to distinguishing features. Similar features do not necessary 
carry identical radiometry values. That uncertainty in either 
manual or automatic interpretation leads to false or missing 
detection. Automatic interpretation methods try to minimize 
the false and missing detection rates. Although a smoothing 
procedure may introduce some false detections, it is also 
necessary to connect some segments of layer features that have 
lost parts and that are still trackable across the echogram. False 
detections, however, need to be minimized. 

D. Example of Application of SAMPA’ s Outputs: 

Accumulation Estimation 

The SAMPA algorithm was applied to a radar echogram 
collected coincidentally with the SEAT 10-1 ice core located 
approximately 20 km from the WAIS divide camp [54]. 
The SAMPA-picked layers were converted from radar travel 
time into fim depth, using densities collected from the ice 
core. The SAMPA algorithm was run over each radar trace 
(« 20 cm along track), and then, peaks in range bin picks were 
averaged to every 100 and 1000 traces to reduce small-scale 
spatial variability, as described in step #5 (see Section IV-E). 
For this application, continuous layers of kilometers are con- 
sidered annual events and not small-scale variability. At core 
site SEAT 10-1, the closest 10 000 traces (« 2 km) were used 
to determine the radar-derived age-depth scale reducing small- 
scale noise. Error bars are ±lcr of the 1000 trace picks, to 
capture the natural layer variation around the core site. The a 
is ±1 standard deviation of the derived accumulation in water 
equivalent (W.E.). Accumulation is derived using the snow 
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SEAT 1 0-1 



Fig. 12. Accumulation rate (in centimeters) of W.E. per year derived from 
radar layers picked using the SAMPA algorithm and derived from the SEAT 10- 
1 ice core [54], showing nearly identical accumulation rates given the estimated 
lcr error. The a is ±1 standard deviation of the derived accumulation in W.E. 
Accumulation is derived using the snow density and radar travel time to the 
SAMPA picks. The mean and standard deviation come from SAMPA picks 
of each radar trace, which are then averaged across the entire radar echogram 
grid, which is « 30 km in length for a single year. The error bars are the same 
distance, i.e., ±lcr in a single year, but the a does change from year to year, 
depending on the SAMPA picks for that year and, ultimately, the snow structure 
(i.e., Some snow layers (or years) have more or less small-scale variability on 
an individual trace level, leading to larger and smaller a values from year to 
year, when averaged over the entire echogram. For instance, a windier year 
would likely have formed sastrugi; therefore, more variability in picks and 
accumulation rates formed point to point (trace to trace) than a calmer year with 
a flatter snow surface that was buried and imaged.) The differences between the 
ice core and the radar are larger on a year-to-year basis than expected. This is 
due to the difference in the type of measurements. 

density and radar travel time to the SAMPA picks. The mean 
and standard deviation come from SAMPA picks of each radar 
trace, which are then averaged across the entire radar echogram 
grid, which is ~ 30 km in length for a single year. The spatial 
average of traces over the echogram ensures that the summer 
season annual layer (dominating layer) is picked and not the 
high- spatial- variability noise from small-scale snow processes. 
The error bars are the same distance, i.e., =b lcr in a single year, 
but the a does change from year to year, depending on the 
SAMPA picks for that year and, ultimately, the snow structure 
(i.e., Some snow layers (or years) have more or less small- 
scale variability on an individual trace level, leading to larger 
and smaller cr values from year to year, when averaged over 
the entire echogram. For instance, a windier year would likely 
have formed sastrugi; therefore, more variability in picks and 
accumulation rates formed point to point (trace to trace) than 
a calmer year with a flatter snow surface that was buried and 
imaged.). The accumulation rate derived from the SAMPA- 
picked layers and the ice core at SEAT 10-1 are nearly identical, 
as shown in Fig. 12. 

The differences between the ice core and the radar are 
larger on a year-to-year basis than expected. This is due to the 
difference in the type of measurements (the radar and ice cores 
are different measurements; the radar is measuring the yearly 
peak in density, whereas the ice core is measuring water isotope 
peak). Further study is needed to understand the difference 
between the yearly peak in density and the isotopic peak, and 
that goes outside the scope of this paper, but over the long term, 
both measurement averages should converge as shown here. 
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VI. Summary 

A SAMPA algorithm for ice- sheet radar echograms has 
been presented. The SAMPA algorithm computes the DRT of 
ice- sheet radar echograms by blocks to overcome undulation 
behavior of firn layers along track. The resulting DRT domain 
carries point sums, i.e., spot peaks representing firn layers, that 
are much easier to extract than in the original space. We applied 
two thresholds (amplitude and width of the peaks) to highlight 
peaks in the RT domain. The backward RT was then computed 
for each corresponding block to get back firn segmented layers 
into their original space. The resulting segmented layers were 
rough and discontinuous. We used a regularization procedure to 
connect and smooth them and to achieve the final layer tracing. 
The overall trackability of good detected layers is greater than 
90%. The SAMPA algorithm has been tested on near-surface 
radar echograms collected by the IceBridge Ku-band radar and 
on ground-based data collected during the 2010 SEAT traverse. 
Firn layers in the area are roughly flat, which leads to making 
use of the RT to extract them. Once the parameters are trained 
and set, SAMPA is fully automated and is capable of tracking 
1000 km of firn layers from radar echograms. 

With respect to the data sets (relatively flat layers) that were 
processed, SAMPA extracted upper layer features better than 
deeper layer features due to low signal returns from attenua- 
tion and subsequent low layers’ contrast below 20 m depth. 
SAMPA-picked layers when compared to the SEAT 10-1 ice 
core show that the algorithm is capable of determining layers 
once trained well enough to give nearly identical accumulation 
rates to a coincident ice core measurement. The SAMPA al- 
gorithm is more efficient for picking multiple fim layers than 
previous manual approaches and can be used for a variety of 
scientific studies, including deriving accumulation rates from 
near- surface radars. 

Our future work will be focusing on reducing false and 
missing detections by 1) managing uncertainty in the detection 
of peaks in the DRT domain via threshold parameters and 
2) by taking into account underlying radiometry values of the 
echogram in the smoothing procedure. 


Appendix I 

Derivation of the Theoretical PDF of the 
Distribution of the RT Domain 

The DRT domain is made with point sum values Si = 
^Sn=i Pn 9 where p n are pixel intensities of the origi- 
nal echogram space, L0 is the length of integration used 
to compute the DRT, and N is the number of elementary 
summations J2n=iPn along a linear feature (see Fig. 1). 
Each p n value follows the log-normal law. C(x\pi,ai) = 
{1 / xcri\/2ii) exp(— ((ln(x) — pi) 2 /2crf)), where pi and cri 
represent the mean and the standard deviation, respectively. 
Therefore, a point sum in the DRT domain follows the sum of 
log-normal distributions. An exact analytical expression of the 
PDF of the sum of log-normal random variables is, however, 
still an unsolved question. 

Using the following assumption: Within a block, pixels be- 
longing to a given segmented layer have close intensity values 
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Po ; therefore, a point sum in the DRT domain can he approx- 
imated by Si = N Yln = 1 Pn ~ N ’ ^0 • po. It becomes much 
easier to assess the PDF h of s*. We have 


Appendix II 

Segmented-Layer Features’ Smoothing 
Algorithm Used in SAMPA 


From the echogram distribution,^ ~ C(jp < 0\pi, (Ji) 
- 1 :C xp 




( 12 ) 


pcriVZn 

From the DRT domain distribution, ~ /i{ Si ,jv,L0,«,<T,}(p) 
= N -L0 £ (iV-LO > O|w ’°0 


1 exp | — n Lx,o) ViY 


paiV2n 


2a? 


03) 


where pi and &i represent the mean and the standard deviation 
of the echogram data, respectively. 

From (5), the probability of detecting peaks Pdp in the DRT 
domain is then given by 


Let si be an extracted segmented-layer feature, which carries 
roughnesses and maybe some missing points due to the absence 
of signal returns or other artifacts. The smoothing algorithm 
[51] estimates the final smoothed layer fl from the extracted 
segmented-layer si using the following iterative scheme: 

fl{k+l} = IDCT(TxDCT(W(sl-fl{k})+fl{k})) (16) 

where 

• DCT and IDCT denote discrete cosinus and inverse 
discrete cosinus transforms, respectively; 

• k represents the current iteration; 

• T a diagonal matrix denoted the tensor defined by 

, if i = j (17) 

0, if * ^ j; 


Uj — 


1 + s (2 — 2 cos(i — l)7r/n) 2 


DRT max 2 

Pdp= J TL(m) 

DRT max i 
DRT max 2 


dm 


~ Adrt — 1 


/ h(m) dm 


DRT inax i 

Ai 


dm (14) 


vdrt 


erf 


^ln( 


DRT m 


N-LO 




- erf 


CTlV 2 y/2J 


<?i 


V2 


V2J 


(15) 


where A^drt is the number of points (i.e., point sum .s, : ) in the 
DRT domain, erf (a;) represents the error function defined by 
erf(x) = (2/ y/n) /* exp(-f 2 ) dt, and [DRT maxl , DRT max2 ] 
is the range of integration. Note that 


" Adrt — 1 



Fig. 3 depicts the theoretical probability of detecting peaks 
Pdp in the DRT domain for Adrt = 50 x 50 point sums, N = 
10, L0 = 5, and the range of point sums containing peaks in 
the DRT domain [DRT max i = 10 9 ,Di?T max 2 = 10 10 ]. On a 
given segmented layer within a block, the Pdp increases as 07 
increases and Pdp > 0.8, for 07 = 10 3 , for the range of point 
sums used for the simulation. This analysis substantiated that 
high signal return echoes (i.e., with high variance of segmented- 
layer points) contribute to a high detectability of peaks in the 
DRT domain. 


• W the bisquare weights defined by [55] 

m = ( (f - ( Mj/4.685) 2 ) 2 , if K/4.685) < 1 
\ 0, if K/4.685| > 1 


where Ui = ri[1.4826MAD(r) 

^-(Vl + A + ^/^Tl+Ife)]- 1 , /•; = *•/,■ ./7 ; 

represents the residual of the Ah observation, MAD 
denotes the median absolute deviation (see [56]), and 

r = rl — fl; 

• GCV ( s ) denotes the generalized cross-validation used to 
estimate the smoothing parameter 5 defined by [57] 


GCV{s ) 


W 1/2 (/f-rQ| 2 /(n-n miss ) 
(1 -Tr(H)/n) 2 


(19) 


where n represents the number of elements of rl ; 

• Tr(H) denotes the matrix trace of H defined by Tr(H) = 
Li[l + s (2 — 2 cos(i — l)7r/n)] -1 , where H = (7 n + 
sD T D)- 1 denotes the hat matrix, and D a tridiagonal 
square matrix defined by [58] as 1 = 2/(/u_i(/u_i + 
T'i))? Di,i — and Di—\{ — ^ j (fiifhi— 1 T 

hi)), hi = \A + y/1 + 16s/(s/2y/l + 16s) denotes the 
step between fl{k} and fl{k + 1}. 

In practice, less than ten iterations are sufficient to get an 
acceptable smoothing result. 
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